Flight planning

In this file, the goal is to quantify the historical cloud cover in the planned flight boxes and to assess the implications for BioSCape operations.

# Load required packages

# Load packages
  library(rgee)
  library(targets)
  library(sf)
  library(tidyverse)
  library(lubridate)
  library(leaflet)
  library(ggplot2)
  library(leafem)



#Load required data

  # get domain
    domain <- st_read("data/output/domain.gpkg",
                      quiet = TRUE)
    domain_sf <- domain
  
  # get flight boxes
    boxes <- st_read("data/flight_planning/20221026_flightboxes.gpkg",
                     quiet = TRUE)
    boxes$id <- 1:20 # need a unique ID to make things easier
    boxes_sf <- boxes
    
  # Download table from drive (to see the code underlying this or to update the data, see the file "R/mock_flights_earth_engine.R")
    googledrive::drive_download(file = "EMMA/cloud_stats.csv",
                                path = "data/test_cloud_stats.csv",
                                overwrite = TRUE) #note: if dling more than one need to add a prefix or something
    
    cloud_table <- read.csv("data/test_cloud_stats.csv")

  # Parse date

    cloud_table %>%
      mutate(year = year(date),
             month = month(date),
             day = day(date),
             day_of_year = yday(date)) -> cloud_table

Flight Boxes

ggplot()+
  geom_sf(data = domain)+
  geom_sf(data = boxes_sf,fill="light green")+
      geom_sf_text(data = boxes_sf,aes(label = id),size=10)+
  labs(x=NULL,y=NULL)
Figure 1. The full domain (light grey) and planned flight boxes (green). Numbers shown are box IDs (used below).

Figure 1. The full domain (light grey) and planned flight boxes (green). Numbers shown are box IDs (used below).

Cloud cover over time

Focal study sites have been selected on the basis of the following criteria:

  #ID by day of year
  
    cloud_table %>%
      group_by(id,day_of_year)%>%
      ggplot() +
      geom_tile(mapping = aes(x = day_of_year,
                              y = id,
                              fill = mean))+
      scale_fill_gradient(low = "sky blue",
                          high = "white")+
      facet_wrap(~year)+
      labs(fill="Mean\nCloud\nCover",
           y="ID",
           x = "Day of Year")
Figure 2. Rows represent different flight boxes (see Fig 1.), columns different days.

Figure 2. Rows represent different flight boxes (see Fig 1.), columns different days.

Mean monthly cloud cover

# mean monthly cloud cover
    cloud_table %>%
      group_by(id, month)%>%
      filter(month != 9)%>%
      summarize(mean_cc = mean(na.omit(mean)))%>%
      inner_join(x = boxes_sf)%>%
      ggplot(mapping = aes(fill = mean_cc))+
      geom_sf()+
      geom_sf(data = domain_sf,
              inherit.aes = FALSE,fill=NA)+
      scale_fill_gradient(low = "sky blue",high = "white")+
      facet_wrap(~month,ncol = 1)+
      geom_sf_text(aes(label = round(mean_cc,digits = 2)))+
      labs(fill="Mean\nCloud\nCover")+
      xlab(NULL)+
      ylab(NULL)
Figure 3. The full domain (light grey) and planned flight boxes, colored by mean monthly cloud cover.Numbers on flight boxes refer to mean cloud cover.

Figure 3. The full domain (light grey) and planned flight boxes, colored by mean monthly cloud cover.Numbers on flight boxes refer to mean cloud cover.

Proportion of clear days for each flight box

Here, we define a “clear” day as one with a mean cloud cover of less than 10 percent.

    #prop clear days
    cloud_table %>%
      na.omit()%>%
      filter(month != 9)%>%
      mutate(binary_clear = dplyr::if_else(mean <= .1,true = 1,false = 0)) %>%
      group_by(id)%>%
      summarize(prop_clear = sum(binary_clear)/n(),
                clear_days = sum(binary_clear),
                total_days = n())%>%
      inner_join(x = boxes_sf)%>%
      ggplot(mapping = aes(fill = prop_clear))+
      geom_sf()+
      geom_sf(data = domain,inherit.aes = FALSE,fill=NA)+
      scale_fill_gradient(low = "white",high = "sky blue",limits=c(0,1))+
      geom_sf_text(aes(label = round(prop_clear,digits = 2)),size=10)+
      labs(fill = "Prop.\nClear",
           x=NULL,
           y=NULL)
Figure 4. Proportion of clear days (mean cloud cover < 10% ) for each flight box. Values in boxes are the proportion of clear days for that box

Figure 4. Proportion of clear days (mean cloud cover < 10% ) for each flight box. Values in boxes are the proportion of clear days for that box

Interactive map

The map below is interactive and allows zooming, moving, and toggling layers on/off.

#Pull in other bioscape layers

    cloud_table %>%
      na.omit()%>%
      filter(month != 9)%>%
      mutate(binary_clear = dplyr::if_else(mean <= .1,true = 1,false = 0)) %>%
      group_by(id)%>%
      summarize(prop_clear = sum(binary_clear)/n(),
                clear_days = sum(binary_clear),
                total_days = n())%>%
      mutate(prop_clear = round(prop_clear,2))%>%
      inner_join(x = boxes_sf)%>%
    st_transform(crs = st_crs(4326)) -> flights

team_requests <- st_read("data/manual_downloads/BIOSCAPE_proposed/20221014_team_polygons.gpkg",quiet = TRUE) %>%
    st_transform(crs = st_crs(4326))

domain_sf %>%
    st_transform(crs = st_crs(4326)) -> domain_wgs84
  

#Make a palette
pal <- colorNumeric(palette = colorRamp(c("white", "blue"), interpolate = "spline"),
                    domain = unique(flights$prop_clear))


#Make labels
labels <- sprintf(as.character(flights$prop_clear)) %>%
  lapply(htmltools::HTML)

leaflet(data = flights) %>%
  addProviderTiles("Esri.NatGeoWorldMap", group = "NatGeo") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "World Imagery") %>%
  addMapPane("flights", zIndex = 420) %>%
  addMapPane("requests", zIndex = 410)%>% 
     addPolygons(stroke = TRUE,
              group = "Flights",
              color = ~pal(prop_clear),
              opacity = 1,
              label = labels,
              options = pathOptions(pane = "flights"))%>%
  addPolygons(data = domain_wgs84,
              stroke = TRUE,
              color = "grey",
              fill = FALSE,
              weight = 3) %>%
      addPolygons(data = team_requests%>%
                st_zm(drop = T, what = "ZM"),
                  stroke = TRUE,
                  color = "brown",
                  group = "Requests",
              options = pathOptions(pane = "requests"))%>%  
    addMouseCoordinates() %>%
    #addImageQuery(sampling_options_wgs84, type="mousemove", layerId = "park_name") %>%
  leaflet::addLegend(position = "bottomright",
            pal = pal,
            values = unique(flights$prop_clear),
            opacity = 1,
            title = "Cluster") %>%
    addLayersControl(
    baseGroups = c("World Imagery","NatGeo"),
    overlayGroups = c("Flights","Requests"),
    options = layersControlOptions(collapsed = FALSE),
    position = "topright")

Figure 5. Interactive Map

Campaign Simulations

      #readRDS("data/temp/sim_output_10pct.RDS")%>%
      readRDS("data/temp/sim_output_05pct.RDS")%>%
        #readRDS("data/temp/sim_output_01pct.RDS")%>%  
        na.omit()%>%
        group_by(start_date)%>%
        summarise(sites_done = sum(!is.na(box_id)),
                  mean_cc = mean(na.omit(cloud_cover)),
                  days_taken = max(date)-min(start_date)+1)%>%
        ungroup()%>%
        mutate(day_of_year = yday(as_date(start_date)),
               year = year(as_date(start_date)))%>%
        group_by(day_of_year)%>%
        summarise(q0 = quantile(days_taken,probs=0),
                  q25 = quantile(days_taken,probs=0.05),
                  q50 = quantile(days_taken,probs=0.5),
                  q75 = quantile(days_taken,probs=0.95),
                  q1 = quantile(days_taken,probs=1)
        ) %>%
        mutate(date = as.Date(day_of_year,origin="2022-12-31"))%>%
        mutate(start_of_month = floor_date(date,unit = "month"),
               month_label = month(start_of_month,label = TRUE),
               julian_label = yday(start_of_month),
               day_of_month = mday(date),
               md_label = paste(month_label,"-",day_of_month,sep = ""))->test
      
      test %>%    
        ggplot()+
        geom_ribbon(aes(ymin=q0,ymax=q1, x = day_of_year),col="grey",alpha=0.2)+
        geom_ribbon(aes(ymin=q25,ymax=q75, x = day_of_year),col="grey",alpha=0.5)+
        geom_line(aes(x=day_of_year,y=q50))+
        geom_hline(yintercept = 37,lty=2)+
        scale_x_continuous(breaks = test$day_of_year,
                           labels = test$md_label)+ #note: it seems like it should be possible to inherit these
        ylab("Flight Days Needed")+
        xlab("Campaign Starting Day")+
        geom_text(label = "estimated total number of flight days",
                  y=37.3,
                  x=280)+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Figure 5. Estimated length of time needed to finish campaign. The dark grey area represents the 50% CI, the light gray the 90% CI, and the solid line represents the median.

Figure 5. Estimated length of time needed to finish campaign. The dark grey area represents the 50% CI, the light gray the 90% CI, and the solid line represents the median.